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Abstract 

We use molecular dynamics (MD) to study the kinetics of surface 
enrichment (SE) in a stable homogeneous mixture (AB), placed in 
contact with a surface which preferentially attracts A. The SE pro- 
files show a characteristic double-exponential behavior with two length 
scales: which rapidly saturates to its equilibrium value, and £ + , 
which diverges as a power-law with time (£+ ~ t e ). We find that hy- 
drodynamic effects result in a crossover of the growth exponent from 
9 ~ 0.5 to 6 ~ 1.0. There is also a corresponding crossover in the 
growth dynamics of the SE-layer thickness. 

1 Introduction 

There is a wide range of physical phenomena associated with binary (AB) 
mixtures in contact with a surface (S). Typically, the surface has a preferen- 
tial attraction for one of the components of the mixture (say, A). A rich phase 
diagram arises, depending on the interactions between AB, AS, BS and the 
temperature T [TJ [2j [31 HJ |5] . For simplicity, we focus on a semi-infinite geom- 
etry. If the bulk system is below the miscibility gap, the AB interface meets 
the surface at a contact angle 9, which is determined by Young's equation 

®- 

"i AB COS# = -Ibs -J AS, (1) 

where 7^/3 denotes the surface tension between a and (3. For jab > 1bs—Jas, 
the surface undergoes partial wetting. For jab < Ibs^Ias, Young's equation 
does not have a solution and the surface is completely wet in equilibrium, 
i.e., the 5-rich phase is expelled from the surface. On the other hand, if the 
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bulk system is above the miscibility gap, there will be no macroscopic phase 
separation. Nevertheless, the surface is enriched in the preferred component 
- the degree of enrichment depends upon various interaction strengths. Both 
the problems of wetting and surface enrichment are of great scientific and 
technological importance. 

We have a long-standing interest in the kinetics of mixtures at surfaces. 
Again, a rich set of problems arises in this context. Consider a homogeneous 
AB mixture (at high temperatures, T = oo) placed in contact with a flat sur- 
face (located at z — 0) which prefers A. If the system is suddenly quenched 
below the miscibility gap at time t — 0, it undergoes phase separation in the 
bulk [HE], and segregates into A-rich and £>-rich domains. Simultaneously, 
the surface is wetted by A. The interplay of these two kinetic processes, i.e., 
phase separation and wetting, results in surface-directed spinodal decom- 
position (SDSD) waves which originate at the surface and propagate into 
the bulk. There has been intense experimental [91 [TOj [HI H2] and theoretical 
[13l[TU[T5] interest in SDSD. Alternatively, we can consider quenches to tem- 
peratures above the miscibility gap or to the metastable region of the phase 
diagram [TB|. In this case, the mixture is stable and continues to be homo- 
geneous for large z. However, the surface becomes enriched in the preferred 
component, resulting in a time-dependent surface-enrichment (SE) profile, 
which propagates into the bulk. This phenomenon of SE is well known in 
the context of polymer blends [I7J [THJ [THl EDI [2H 122] and biopolymer mixtures 
[23] > binary alloys [2U [2S] , and the wetting of surfaces of fluids [U |2S] • 

In this paper, we undertake a comprehensive molecular dynamics (MD) 
study of the kinetics of SE. This paper has two primary goals. First, we would 
like to test the theoretical results for enrichment kinetics (summarized later) 
via a realistic microscopic MD simulation. The available numerical studies 
of this problem have focused on Langevin simulations of phenomenological 
models. Second, we would like to examine whether hydrodynamic effects 
(which are naturally incorporated in our MD simulations) have any effect on 
SE kinetics. After all, many experiments in this area involve fluid mixtures or 
polymer mixtures, where velocity fields play a significant role. In the context 
of phase-separation kinetics, we know that hydrodynamic effects drastically 
alter the late-stage dynamics - both without surfaces [H [7J, [23 [23 129] and 
with surfaces [H EU E2] ■ 

This paper is organized as follows. Section [2] describes the details of our 
MD simulations. In Sec. [3J we present comprehensive numerical results and 
compare them with theoretical predictions. Finally, Sec. H] concludes this 
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paper with a summary and discussion. 



2 Details of Molecular Dynamics Simulations 

Following Das et al. [33], we consider a binary (AB) fluid mixture of point 
particles confined in a rectangular box of volume V = L x L x D. We apply 
periodic boundary conditions in the x and y directions, while impenetrable 
walls or surfaces are present at z = and z = D. These surfaces give rise to 
an integrated Lennard- Jones (LJ) potential (a = A,B): 
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where n is the reference density of the fluid, and a is the LJ diameter of the 
particles. In Eq. e r and e a are the strengths of the repulsive and attractive 
parts of the wall potentials. For the wall at z = 0, we choose 8^ = 1 and 
5b = 0, so there is only a repulsion for B particles, while A particles are 
attracted at large distances and repelled at short distances. For the wall 
at z = D, we set 5 a = 5 b — 0, so that both A and B particles feel only 
repulsion. Further, z! = z + a/2 for the wall at z — 0, and z' = D + a/2 — z 
for the wall at z — D. Therefore, the singularities of u w (z) do not occur 
within the box range < z < D, but rather at z — —a/2 and z = D + a/2, 
respectively. 

The particles in the system interact with LJ potentials: 



u(rij) = 4e Q/3 



a 
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where = |rj — f^|; a, (3 = A, B. The energy scales are 



(3) 
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for which the equilibrium phase behavior is well studied [MJ |35j [36] . We use 
the truncated, shifted and force-corrected LJ potential with = r c = 2.5a 
[37]. We consider the case with equal numbers of A and B particles (Na = 
N B = N/2), and their masses are set to be equal, rriA — tub — m — 1. We 
also set a — 1, e = 1 and ks — 1, such that the MD time unit is 
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We work with a high-density liquid having n = N/V = 1, which makes 
the system incompressible. Notice that crystallization is not a problem at 
the temperatures of interest. Finally, simulations were run for three surface 
energy strengths - e a — 1.2, 2.0, 3.0 with e r = 0.5. 

For our study, we chose L = 32 and D = 64 (N = 66536 particles). 
As the bulk remains homogeneous, the lateral size L (in the x, y directions) 
is not severely constrained. However, in the direction perpendicular to the 
surface located at z — 0, we need sufficiently large D to ensure decay of the 
enrichment profiles as z — > D. For the range of times studied here (t < 7000), 
test runs with other linear dimensions showed that D = 64 is large enough to 
eliminate finite-size effects, within the limits of our statistical accuracy. For 
a smaller system size (D = 32), we encounter finite-size effects by t ~ 3000. 
The statistical quantities presented here were obtained as averages over 50 
independent runs. The MD runs were carried out using the standard Verlet 
velocity algorithm [38J, with a time-step At = 0.07 in MD units given by 
Eq. ([5]), and the Nose-Hoover thermostat [38]. We prepared the homogeneous 
initial state for a run by equilibrating the mixture of N particles at high T, 
in the specified volume with periodic boundary conditions in all directions. 
At time t = 0, the system is quenched to T = 2.0 ~ 1.41T C (T c ~ 1.423) 
[35l 136] . and surfaces are introduced at z = 0, D. 

3 Theoretical Background and Numerical Re- 
sults 

We are interested in the time-dependent morphology which arises during SE. 
We characterize the morphology via laterally-averaged depth profiles and 
their various properties, e.g., surface value of order parameter, zero-crossings, 
moments, etc. Before proceeding, it is useful to summarize theoretical results 
in this context. 

3.1 Theoretical Background 

Jones et al. [HI US] have studied the kinetics of SE for polymer mixtures, and 
found that the enrichment profiles are characterized by diffusive length scales. 
Numerical studies of this problem (with both short-ranged and long-ranged 
surface fields) have been conducted by Jiang and Ebner [39] using MC simu- 
lations, and Toral and Chakrabarti [10] via Langevin simulations. Binder and 
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Frisch [JI] studied SE in a linearized Ginzburg-Landau (GL) model with a 
delta-function (short-ranged) surface potential. Their GL model was derived 
by coarse-graining the Kawasaki spin-exchange kinetic Ising model in the 
presence of a surface. Their results are expected to be valid in the weak-field 
regime, where the concentration variations are sufficiently weak that nonlin- 
ear effects can be neglected. Puri and Frisch (PF) [42J have also studied the 
case with a short-ranged surface potential in the framework of a linear theory. 
They also undertook Langevin simulations of the corresponding GL model 
for both weak and strong surface fields. In the latter case, linear theory is 
no longer applicable as the order parameter value at the surface becomes ap- 
preciably large. PF demonstrated that the morphological features predicted 
by linear theory are also seen in the nonlinear regime. These results were 
extended to the case with an arbitrary surface potential V(z) by Frisch et 
al. (FPN) [43] . Here, we briefly summarize the results of PF and FPN. 

Let us first discuss the modeling of segregation kinetics at surfaces [HI [15]. 
The bulk order parameter satisfies the Cahn-Hilliard (CH) equation: 







1. 



-0(p,z,t) = -V 2 [sgn(T c -T)0-0^ + -V 2 0-^)j, z>0. (6) 

Here, all the quantities have been rescaled into dimensionless units [B] . The 
order parameter 0(r, t) ~ n^r, t) — riB(r,t), where n a (r, t) denotes the den- 
sity of species a at space-point r and time t. We have decomposed coordinates 
as r = (p,z), where p and z denote coordinates parallel and perpendicular 
to the surface located at z = 0. The function sgn(x) = 1 for x > and — 1 
for x < 0. The surface potential V(z) (< 0) is chosen so that it enriches the 
surface in A. 

Equation (Q must be supplemented by two boundary conditions at z — 0, 
as it is a fourth-order partial differential equation. Since the surface value of 
the order parameter is not a conserved quantity, we assume a nonconserved 
relaxational kinetics with time-scale r for this quantity at the surface: 



d d 
To 4>{P, 0, t) = h 1 + g<p(p, 0, t) + 7^0(P, *, *) 



+ 7Vy0(p,O,t), (7) 

2 = 



where h% = — V(0), and (7,7,7 are phenomenological parameters, which are 
related to the bulk correlation length [H] . Finally, we implement the no-flux 
boundary condition at the surface, which enforces order parameter conserva- 
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tion: 







d_ 

dz 



sgn(T c - T)0 - 3 + lv 2 - V(z) 



z=0 



(8) 



Equations ([H])-® describe both SDSD (when T < T c ) and SE (when T > T c ), 
as long as the dynamics is diffusive. This is appropriate for phase separa- 
tion in solid mixtures or the early stages of segregation in polymer mixtures. 
However, most experiments involve fluid mixtures, where hydrodynamic ef- 
fects play an important role in the intermediate and late stages of phase 
separation. At a phenomenological level, hydrodynamic effects can be incor- 
porated via the Navier-Stokes equation for the velocity field - the resultant 
coupled equations are known as Model H [201 US]- This must be supple- 
mented by appropriate boundary conditions at the surfaces. Alternatively, 
one can consider microscopic models of fluid mixtures at a surface, which 
naturally incorporate the fluid velocity field. We adopt the latter strategy in 
this paper, and use MD simulations to study SE in fluid mixtures. 

For reference, it is useful to summarize results for the diffusive case. We 
consider the model in Eqs. ©-flB]) for T > T c , as we are interested in the 
kinetics of SE. In this case, the order parameter field remains homogeneous 
in the direction parallel to the surface, i.e., we can neglect the p -dependence 
of the order parameter <p(p,z,t) ~ <p(z,t). FPN have solved the linearized 
version of this model, which is appropriate when is small and the 3 - 
term can be neglected in Eqs. fl6]) and flH]). Seeing that the bulk remains 
homogeneous, there is a significant enhancement of <fi only for z ~ 0. Thus, 
the linearized model is valid for weak surface fields. In this case, the results 
of FPN are as follows. 

The SE profiles have a double-exponential form |41j : 



where the amplitudes B_(t), B + (t) > 0. Notice that the conservation con- 
straint dictates that -B_£_ = B + £ + . The quantities B_(t) and rapidly 
saturate to their equilibrium values which depend on the surface potential: 



The other length scale grows diffusively with time, and B + (t) shows a 



(f)(z,t) ~ B_(t) e-*K-® - B+(t) e^+W, 



(9) 



(10) 
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corresponding decay: 



B+(t) ~ ait' 1 ' 2 , 

£f(f) ^ M 1/2 - (11) 

The other properties of the enrichment profiles are obtained from Eq. ([9]). 
An experimentally important quantity is the time-dependence of <\>{z = 0,t), 
the value of the order parameter at the surface [201 Ell E2] : 

0(0, £) ~Oi -a 2 t~ 1/2 . (12) 

Thus, 0(0, £) saturates diffusively to its equilibrium value. The thickness of 
the enrichment layer is measured as the first zero of the double-exponential 
profile. This increases logarithmically with time: 

- H9- (l3) 

Finally, consider the time- dependence of the profile moments: 

/•OO 1 

(z m ) = dzz m <P(z,t) ~ — - ) • (14) 

Therefore, the asymptotic (i — >■ oo) behavior of (z m ) is (z m ) ~ t m l 2 . 

These results apply for diffusive transport and are universal for a wide set 
of potentials. They have been obtained in the context of a linear theory, and 
there is a range of weak surface fields where the analytical results obtained 
from the linear model agree well with the numerical solution of the nonlinear 
Eqs. (jBJ)-(jSi) 03]- As the field strength is increased, the validity of linear 
theory breaks down in the vicinity of the surface as the degree of enrichment 
becomes larger. Nevertheless, FPN demonstrated that the diffusive behavior 
of various profile characteristics is unaffected, even in the strongly nonlinear 
regime. 

As stated earlier, we have undertaken MD simulations to examine whether 
hydrodynamic effects have any impact on the above phenomenology. Let us 
next present results from these simulations. 



7 



3.2 Numerical Results 



Figure [T] shows a three-dimensional snapshot of SE in a binary (AB) mixture 
at t — 7000. The surface field strengths are e a = 3.0 and e r = 0.5 in 
Eq. 02]). We see the formation of an A-rich (marked in gray) layer at the 
surface (z = 0), resulting in a time-dependent SE profile which propagates 
into the bulk. However, in the bulk (large z), the thermodynamically stable 
mixture continues to be homogeneous. In Fig. |2j we show cross-sections of 
the snapshots at z = for t = 70, 7000. The cross-section shows all A atoms 
(marked gray) and all B atoms (marked black) lying in the interval z G [0, a]. 

In Fig. El we show the temporal evolution of the laterally-averaged order 
parameter profiles av (2;,t) vs. z, obtained from our MD simulations [44J. 
The order parameter is defined in terms of the local densities as 0(r, t) = 
(jiA —tib)/ (ua + tib)- The quantity av (^, t) is obtained by averaging 0(r, t) 
in the directions parallel to the surface, and then further averaging over 50 
independent runs. These laterally-averaged profiles are analogous to depth 
profiles measured in experiments - see Fig. 7 in Ref. [20] or Fig. 4 in Ref. [22J. 
The enrichment profiles are shown for the case with field strength e a = 2.0 at 
times t = 280, 1400, 7000. It is clear that a layer rich in A forms at the surface 
immediately after the field is turned on. Due to the conservation of the order 
parameter, there must be a corresponding depletion layer which decays to 
av ~ in the bulk. These profiles are in agreement with the experimental 
observations of Jones et al. [TS] on blends of deuterated and protonated 
polystyrene, and the experimental results of Mouritsen [23J on biopolymer 
mixtures. Notice that similar profiles are seen for SDSD or surface-directed 
phase separation if the system is quenched to the metastable region of the 
phase diagram [16J. The evolution dynamics in that case is analogous to the 
SE problem as long as droplets are not nucleated in the system. 

We make the following observations about the enrichment profiles in 
Fig. |3j First, for the field strengths we consider, the surface is strongly 
enriched in A: <p av (0,t) ~ 1 for e a = 2.0,3.0. Thus, the linear theory is not 
applicable in the enrichment layer as this would require 3 <C 0. Second, 
at late times (say t = 7000 in Fig. [3]), the depletion region stretches deep 
into the bulk. To avoid finite-size effects, we confine our simulation to time 
regimes where the thickness of the depletion region CD = 64, the box size 
in the ^-direction. As shown in Fig. [31 the profiles are fitted very well by 
the superposition of two exponential functions as in Eq. fl9]). Our simula- 
tions show that B_{t) and £_(t) rapidly saturate to their equilibrium values 
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(in agreement with the predictions of linear theory j?U H5J H3]) and may be 
treated as static parameters. 

In Fig. HI we show the time-dependence of £_ (t) and (t) for three surface 
field strengths, e a = 1.2,2.0,3.0. The saturation value of £_ increases with 
the field strength. Further, we find that grows with time as a power- 

law (£ + ~ t 9 ) but there is a clear crossover in the growth exponent. At early 
times (t <C t c ), we have 9 ~ 0.5, in the conformity with the prediction of 
linear diffusive growth. However, there is much more rapid growth at late 
times (t ^> t c ) with 9 ~ 1.0. To understand the crossover in Fig. HJ consider 
the dimensionless evolution equation for the order parameter in the presence 
of a velocity field v(f, t) [TJ |S]: 

d 

—<j>{r,t) = V 2 fi-v-V(p, (15) 



dt 



V = + 3 -^V 2 < 



For the SE problem, we set 4>{f,t) ~ <f>(z,t) and v z (f, t) ~ v z (z,t). Then 

We use the double-exponential form of (j>(z,t) in Eq. (j9]) to estimate the 
various terms in Eq. ffl6l) to leading order at z ~ 0(^ + ), i.e., far from the 
surface. We have 

90 1 ^+ 



0t ^ rft 

i 

r\j 

dz 2 a ' 



dz 



-2 • 



(17) 



We have used the general relation = to obtain the above ex- 

pressions. In this case, the bulk is homogeneous and there is no structure 
formation in the density or velocity fields, i.e., v z ~ constant. At early times, 
the diffusive term in Eq. ffl6|) dominates, yielding £ + ~ t 1 ! 2 . At late times, 
the convective term in Eq. f fl6|) is dominant, giving £ + ~ v z t. The precise 
dependence of the crossover on various physical parameters can be estimated 
by considering the dimensional version of Eq. (|15|). 
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The crossover in the growth exponent 9 is reminiscent of phase-separation 
kinetics in fluid mixtures quenched below T c [HE]- I n that case, the domains 
grow as L(t) ~ t x with the exponent crossing over from 1/3 (diffusive regime) 
to 1 (viscous hydrodynamic regime) to 2/3 (inertial hydrodynamic regime). 
Recently, we have observed the 1/3 — >■ 1 crossover in the growth law for 
the wetting layer in MD studies of surface-directed phase separation [32J. 
Our MD results in this paper show that convective transport accelerates the 
growth of the SE layer also. Further, the onset of the hydrodynamic regime 
is faster for stronger surface fields. This result has important experimental 
implications, and we urge experimentalists to undertake a detailed study of 
this problem. 

We have also confirmed that B_(t) rapidly saturates to its equilibrium 
value, and that B + {t) ~ £ + (£) _1 (results not shown here). The profile pa- 
rameters are consistent with the conservation constraint, -B_£_ — B + ^ + . 
Figure |5] shows the time- dependence of the surface value of the order param- 
eter. We plot av (O, oo) — av (O, t) vs. t^ 1 for e a = 1.2, 2.0, 3.0, demonstrating 
that 0av(O, t) saturates linearly to its asymptotic value av (O, oo). Recall that 
0(0, t) ~ B_{t) — B + {t) ~ ai — a 2 /£+, and ~ t in the late stages. 

Finally, we focus on the time-dependence of the thickness of the enriched 
layer, which is also of considerable experimental interest. For the double- 
exponential profile in Eq. fl9]), the zero is located at [cf. Eq. (TT3]) ] 



Thus, we expect Z (t) ~ In t in both the time- regimes of Fig. HI but the slope 
should be steeper for t > t c . This is precisely the behavior seen in Fig. |6l 
where we plot our MD results for Zo(t) vs. t on a log-linear scale. This 
confirms that there is a crossover in the growth exponent from the diffusive 
regime (8 ~ 0.5) to the hydrodynamic regime (6 ~ 1.0). 

4 Summary and Discussion 

Let us conclude this paper with a brief summary and discussion of our results. 
We are interested in the kinetics of surface enrichment (SE), which occurs 
when a miscible or metastable binary (AB) mixture is placed in contact with 
a surface having a preferential attraction for one of the components (say, A). 
A closely-related problem is that of surface-directed spinodal decomposition 




(18) 
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(SDSD), where an unstable homogeneous mixture is placed in contact with 
a wetting surface [2]. The problems of SE and SDSD are of great scientific 
and technological importance. Experiments in this area have been performed 
on polymer blends, fluid mixtures, alloys, etc. 

In this paper, we undertake comprehensive molecular dynamics (MD) 
simulations to study the kinetics of SE. The typical enrichment profile con- 
sists of an enriched surface layer, followed by an extended shallow depletion 
region. This profile propagates into the bulk with the passage of time. We 
are interested in understanding the role of hydrodynamic effects in driving 
the growth of the enrichment layer. In the context of phase-separation kinet- 
ics, we know that hydrodynamic transport drastically alters the intermediate 
and late stages of domain growth. At early times, the characteristic length 
scale of the SE profile grows diffusively with time (£+ ~ t 1 ^ 2 ), which is con- 
sistent with linear theory for the Cahn-Hilliard model. However, the growth 
exponent undergoes a crossover to a convective regime, and the late-stage 
dynamics is £ + ~ t. There is a corresponding crossover in the growth dy- 
namics of the thickness of the enrichment layer Z (t). The growth of Z (t) 
is logarithmic in both regimes, but with a different slope. 

The MD results presented here have significant implications for SE ex- 
periments, as many of these are performed on fluid mixtures. We hope that 
our MD results will provoke fresh experimental interest in this problem, and 
our theoretical results will be subjected to an experimental confirmation. 
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t = 7000 




Figure 1: Snapshot of surface enrichment (SE) in a binary (AB) Lennard- 
Jones (LJ) mixture, which is confined in a box of size Lx Lx D with L = 32, 
D = 64. (All lengths are measured in units of the LJ diameter.) Periodic 
boundary conditions are applied in the x, y directions, while an impenetrable 
L x L surface at z = attracts the A-particles. For clarity, we show only 
part of the simulation box with z G [0,32]. The initial condition for this 
run consisted of a random mixture of equal amounts of A and B particles 
(Na = Nb = 32768). Time is measured in dimensionless LJ units. For 
further details of the MD simulation, see Sec. [2j The A-particles are marked 
gray, and the -B-particles are marked black. The snapshot corresponds to 
t = 7000 for the surface potential in Eq. (J2]) with e a = 3.0, e r = 0.5. The A- 
rich enrichment layer is formed at z — 0, while for large z, the stable mixture 
continues to be homogeneous. 
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t = 70 



t = 7000 




Figure 2: Cross-section at z = of the SE snapshots at time t = 70 and 7000. 
The A-particles are marked gray and the S-particles are marked black. Other 
simulation details are the same as for Fig. [TJ 
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Figure 3: Laterally-averaged order parameter profiles [(f>g, v (z, t) vs. z] at the 
dimensionless times t = 280, 1400, 7000. The surface potential is given by 
Eq. (J2]) with e a = 2.0, e r = 0.5. The double-exponential fits for the SE 
profiles are shown as solid lines. 
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Figure 4: Time-dependence of the fit parameters £_ (left panel) and £ + 
(right panel) of SE profiles. The surface field values are indicated in the 
legends. The length scale £_(£) rapidly saturates to its equilibrium value, 
whereas £+(£) shows a crossover from the diffusive regime (£+ ~ t 1 / 2 ) to 
the hydrodynamic regime (£+ ~ t) for higher surface field strengths (e a = 
2.0, 3.0). The lines of slope 1/2 and 1 are provided as a guide to the eye. 
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Figure 5: Time-dependence of the surface value of the order parameter for 
the SE profiles. We plot av (O, oo) - <p av (0,t) vs. r 1 for e a = 1.2,2.0,3.0. 
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Figure 6: Time-dependence of the first zero crossing of the SE profiles. We 
plot Zo(t) vs. t (note the logarithmic scale of time) for e a = 1.2, 2.0, 3.0. This 
plot confirms the crossover in the growth exponent from the diffusive regime 
to the hydrodynamic regime. 
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